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In the present work we investigate the main differences in the lead neutron skin thickness, binding 
energy, surface energy and density profiles obtained with two different density dependent hadron 
models. Our results are calculated within the Thomas-Fermi approximation with two different 
numerical prescriptions and compared with results obtained with a common parametrization of the 
non-linear Walecka model. The neutron skin thickness is a reflex of the equation of state properties. 
Hence, a direct correlation between the neutron skin thickness and the slope of the symmetry 
energy is found. We show that within the present approximations the asymmetry parameter for low 
momentum transfer polarized electron scattering is not sensitive to the model differences. 
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I. INTRODUCTION 

The relation between neutron star properties which are 
obtained from adequate equations of state (EoS) and the 
neutron skin thickness has long been a topic of investi- 
gation in the literature. The details of this relation and 
the important quantities to be discussed have been well 
established in [l| , where it was shown that the difference 
between the neutron and the proton radii, the neutron 
skin thickness, is linearly correlated with the pressure of 
neutron matter at sub- nuclear densities. This is so be- 
cause the properties of neutron stars are obtained from 
appropriate EoS whose symmetry energy depends on the 
density and also controls the size of the neutron skin 
thickness in heavy and asymmetric nuclei, as 208 Pb, for 
instance. It is important to remember that the EoS in 
neutron stars is also very isospin asymmetric due to the 
/3- equilibrium constraint. 

Hence, isospin asymmetry plays a major role in the un- 
derstanding of the density dependence of the symmetry 
energy and the consequences it may arise 0- In [3, 0] 
it was shown that the models that yield smaller neutron 
skins in heavy nuclei tend to yield smaller neutron star 
radii due to a softer EoS. 

Neutron stars are believed to have a solid crust formed 
by nonuniform neutron rich matter in /3-equilibrium 
above a liquid mantle. In the inner crust nuclei coex- 
ist with a gas of neutrons which have dripped out. The 
properties of this crust as, for instance, its thickness and 
pressure at the crust-core interface depend a lot on the 
density dependence of the EoS used to describe it 0, @ . 
On the other hand, it is well known [|| 0] that the ex- 
istence of phase transitions from liquid to gas phases in 
asymmetric nuclear matter (ANM) is intrinsically related 
with the instability regions which are limited by the spin- 
odals. Instabilities in ANM described within relativistic 
mean field hadron models, both with constant and den- 



sity dependent couplings at zero and finite temperatures 
have already been investigated [7( and it was shown that 
the main differences occur at finite temperature and large 
isospin asymmetry close to the boundary of the insta- 
bility regions. In neutral neutron-proton-electron (npe) 
matter the electrons are also included. In a thermody- 
namical calculation the instabilities almost completely 
disappear due to the high electron Fermi energy [8J. 

However, in a dynamical calculation which includes the 
Coulomb interaction and allows for independent neutron, 
proton and electron fluctuations 0, [l(| , it is seen that 
the electron dynamics tends to restore the short wave- 
length instabilities although moderated by the high elec- 
tron Fermi energy. 

Moreover, it is also known that the liquid-gas phase 
transition in ANM can lead to an isospin distillation phe- 
nomenon, characterized by a larger proton fraction in 
the liquid phase than in the gas phase. This is due to 
the repulsive isovector channel of the nuclear interaction 

In a recent work the spinodal section and related quan- 
tities, as the neutron to proton density fluctuations re- 
sponsible for the distillation effect, has been studied 
within different relativistic models [8| . It was shown that 
the distillation effect within density dependent relativis- 
tic models decreases with density above a nuclear density 
of ~ 0.02 — 0.03 fm -3 , a result similar to the one obtained 
with the SLy230a parametrization of Skyrme interaction 
[l4| and contrary to the results found with the more com- 
mon relativistic parametrizations with no density depen- 
dent coupling parameters. In the last case the distillation 
effect becomes always larger as the density increases. 

Also, the behavior of the symmetry energy obtained 
with density dependent models is closer to what one ob- 
tains with non-relativistic models than with other rel- 
ativistic models with constant couplings p}. In an at- 
tempt to understand this behavior, a comparison be- 
tween the non-relativistic Skyrme effective force and rel- 
ativistic mean field models at subsaturation densities was 
performed jl5| . It was shown that the relativistic mod- 
els could also be reduced to an energy density functional 
similar to the one describing the Skyrme interaction. 
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There have already been some efforts in order to com- 
pare nuclear matter and finite nuclei properties obtained 
both with relativistic and non-relativistic models 

HED3 

but there is no clear or obvious explanations for the differ- 
ences. At very low densities both, the relativistic and the 
non-relativistic approaches predict a non-homogeneous 
phase commonly named pasta phase, formed by a com- 
petition between the long-range Coulomb repulsion and 
the short-range nuclear attraction [l8l ]. 

Based on the above arguments, it is very important 
that an accurate experimental measurement of the neu- 
tron skin thickness is achieved. This depends on a precise 
measurement of both the charge and the neutron radius. 
The charge radius is already known within a precision of 
one percent for most stable nuclei, using the well-known 
single-arm and non-polarized elastic electron scattering 
technique as well as the spectroscopy of muonic atoms 
[l9| . For the neutron radius, our present knowledge has 
an uncertainty of about 0.2 fm [201 ] . However, using po- 
larized electron beams it is possible to obtain the neutron 
distribution in nuclei in a fairly model independent way, 
as first discussed in [21[ and, as a consequence, to obtain 
the desired neutron radius. In fact, the Parity Radius 
Experiment (PREX) at the Jefferson Laboratory [22j is 
currently running to measure the 208 Pb neutron radius 
with an accuracy of less than 0.05 fm, using polarized 
electron scattering. 

In the present work, we use two different hadronic 
models that incorporate density dependence in differ- 
ent ways. The first one, to which we refer next as the 
TW model is a density dependent hadronic model with 
the meson-to-nucleon couplings explicitly dependent of 
the density [23|, H3|. In the following it is used to cal- 
culate the neutron skin thickness of 208 Pb, which is a 
neutron-rich heavy nucleus. This model was chosen be- 
cause it is based on a microscopic calculation, fits well 
many nuclei properties and, as stated above, has shown 
to provide results which are different from the usual NL3 
[25| and TM1 [26| parametrizations for the non-linear 
Walecka model (NLWM), having a richer density depen- 
dence of the symmetry energy than most of the rela- 
tivistic nuclear models. The original motivation for the 
development of this density dependent hadronic model 
[27l [28[ was to reproduce results obtained with the rel- 
ativistic Dirac-Brueckner Hartree-Fock (DBHF) theory 
[29l ] . Later the DBHF calculations for nuclear matter 
were taken only as a guide for a suitable parametrization 
of the den sity dep endence of the meson-nucleon coupling 
operators [2J,[30|. Moreover, density dependent hadronic 
models can also be a useful tool in obtaining EoS for 
neutron stars even if hyperons are to be considered [12] , 
which is not the case if NL3 or TM1 are used. Both, NL3 
and TM1, can only be used if the EoS is restricted to ac- 
commodate neutrons, protons and the leptons necessary 
to enforce /3-stability. Once hyperons are included, the 
nucleons acquire a negative effective mass above ~ 3— 4po 
densities [33l . |34| , where po is the nuclear saturation den- 
sity. 



The second model, that we refer to as NLuw model, in- 
cludes non- linear a — p and ui — p couplings [3, [HI, [3(| [3?| 
which allow to change the density dependence of the sym- 
metry energy of the most common parametrizations of 
the NLWM that show essentially a linear behavior of the 
symmetry energy with density. However, the symme- 
try energy determines the behavior of isospin asymmetric 
matter and therefore is intrinsically related to the char- 
acteristics of the EoS that can describe neutron stars. 
Within this model the authors of [|| have shown that 
the neutron skin thickness of 208 Pb was sensitive to the 
isovector channel of the nuclear interaction and there was 
a correlation between neutron skin thickness of nuclei and 
properties of neutron stars. 

For the sake of completeness, the results of the present 
work, whenever possible are compared with the results 
obtained with the NL3 parametrization of the NLWM, 
known to describe finite nuclei properties well. 

We perform two different numerical calculations to ob- 
tain the 208 Pb properties: a Thomas-Fermi approxima- 
tion based on the liquid-gas phase transition developed 
in [3a ] and a Thomas-Fermi approximation based on a 
method proposed in [39] , where a harmonic oscillator ba- 
sis is used. We restrict ourselves to the Thomas-Fermi 
approximation because, as we show in the Results sec- 
tion at the end of the paper, for the purpose of obtaining 
correct surface energy and neutron-skin thickness, it is 
almost as good as the solution of the Dirac equation. 

At this point it is worth mentioning that the scalar- 
isovector S mesons, which play an important role in the 
isospin channel, could also be incorporated in our work 
as done in 0, H, S3] but in order to make the compar- 
isons among different approximations as simple as possi- 
ble, they will be included in a future work. Finally, as 
we are interested in nuclei ground state properties, all 
calculations are performed at zero temperature. 



II. THE TW DENSITY DEPENDENT 
HADRONIC MODEL 

Next we describe the main quantities of the TW model, 
which has density dependent coupling parameters. The 
Lagrangian density reads: 



C = tp 



7m id" - T V V» 



\m 2 W - \b^-B»" + \m 2 p h^ - \F„ V F^ (1) 
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where </>, V^, b M and are the scalar-isoscalar, vector- 
isoscalar and vector-isovector meson fields and the pho- 



ton field respectively, = d^Vu — d v Vp 



B 



<9„b p - r p (b M x b„), = d^Ay - d u A p and 



T p3 = 1; and t„3 = — 1. The parameters of the model 
are: the nucleon mass M — 939 MeV, the masses of the 
mesons m x , m„, m p , the electromagnetic coupling con- 



stant e = y47r/137 and the density dependent coupling 
constants T s , T v and r p , which are adjusted in order 
to reproduce some of the nuclear matter bulk properties 
shown in Table [H using the following parametrization: 



Ti(p) = Ti(p sat )hi(x), x = p/p S at, 



(2) 



with 



and 



where E p ^ is the usual vector self-energy, pu p — j p with 
u 2 = 1 j u = ipjvip, J3 = ipr^"tp and, as a result of 
the derivative of the Lagrangian with respect to p a new 
term appears, £ p , which is called rearrangement self- 
energy and has been shown to play an essential role in 
the applications of the theory. This term guarantees the 
thermodynamical consistency and the energy-momentum 
conservation. For more detailed calculations, at zero and 



finite temperatures, please refer to [41 1. 

In the static case there are no currents in the nucleus 
and the spatial vector components are zero. Therefore, 
the mesonic equations of motion become: 



m\(f> - T s p s , 
■ mlV - T v p, 

V 2 b =m 2 p b -^p^ 



Wo 



(12) 
(13) 

(14) 



hp(x) = exp[-a p (a; - 1)], 



(4) 



with the values of the parameters m^, Fj(p sat ), a^, bi, Ci 
and di, i = s,v,p given in [24!] . This model does not 
include self-interaction terms for the meson fields (i.e. 
k = 0, A = and £ = ) as in NL3 or TM1 parametriza- 
tions for the NLWM. 

The field equations of motion follow from the Euler- 
Lagrange equations. When they are obtained, some care 
has to be taken since the coupling operators depend on 
the baryon fields V> and ip through the density. When the 
partial derivatives of C are performed relatively to the 
fields tp and ip, they yield extra terms due to the func- 
tional dependence of the coupling operators. The new 
terms are absent in the usual Quantum Hadrodynamic 
(QHD, NLWM) models [H, US SH]. The equations of 
motion for the fields read: 





+ ml)(f> 


= r s ^v> 


(5) 


<w 




= r„^yv> 


(6) 


d v B ^ 


+ m 2 p b» 




(7) 






= ^(l + Tah^V, 


(8) 




- M*] ip 


= , 


(9) 



where M* = M—T s 4>. Notice that in the equation of mo- 
tion for the baryon field tp the vector self-energy consists 



of two terms, S, 



_ v (°) 



S^, where: 

H 1 



= r w V^ + • b M + |(1 + r 3 )A„ (10) 



(11) 



-ep P , 



(15) 



where p s —< ipip > is the scalar density, p = p p + p n , 
P3 = Pp — Pn and p P and p n are the proton and neutron 
densities. 



A. Thomas-Fermi approximation 



We first define the functional 

^ = E — PpBp — p n B n , 



(16) 



where E is the energy, p p (p n ) is the proton (neutron) 
chemical potential and B p (B n ) is the proton (neutron) 
number. Within the semi-classical Thomas-Fermi ap- 
proximation, the energy of the nuclear system with par- 
ticles described by the one-body phase-space distribution 
function /(r, p, t) at position r, instant t with momentum 
p is given by 



E =J2' 



d 3 P 



d 3 r [(V</>) 2 



/i(r,p,t) 



M 



2 ±2 
m s <p 



(W ) 2 - m 2 v V 2 



-(Vb ) 2 -m 2 b 2 ~(yA ) 2 



(17) 



where 



v„ = r,i;, + ^6 ,.i„. v„ =r,i;,-if^ 
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7 = 2 refers to the spin multiplicity and the distribu- 
tion functions for protons and neutrons are 



fi = 0(k 2 Fi (r)-p 2 ), 
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In this approach, the scalar, proton and neutron densities 
become: 



Ps(r) 



,fc Fi (r) M * 



2tt 2 



III. NLwp MODEL 

The Lagrangian density that incorporates the extra 
non-linear a — p and u> — p couplings [H, [H, [H, [13 reads 



d 3 rp l; Pl {r) = g-2fc«(r). 



with e = Vp 2 + -^* 2 an d 
B 

From the above expressions we get for ([1 

n = 

with 



£ = V 



7a< 



- e (1 + Tl3) ^ I — (M — g s 



1> 



d 3 r Q [(V0) 2 - (W ) 2 - (V6 ) 2 - (VAq) 2 ] + V e ^j +1(0^8^ - m 2 J 2 ) - - ±\cf> 4 



V ef = - [m 2 s cj) 2 - m 2 v V£ - m 2 p bl] - p p p p - p n p n 



2tt 2 



i—p,n 



I dpp 2 e + r v V p + T p ^-p 3 + eA p p (18) 



Minimization of f2 with respect to kpi(r), i = p,n, 
gives rise to the following conditions 



k l P - ^k 2 Fp + M* 2 - T V V - ^-b - eA Q - ) = 
and 



k 2 Fn Un - V4„ + M* 2 - FvVo + -fb Q - Efj = 0, 
where the rearrangement term is 

R dv v dv p b or s 

X =-Q i PV + -g i P 3 ^--g i P S t. 

From the above equations we obtain kp p — and kp n = 
or, for kF P or kp n different from zero, 



Pp = \/ k2 F P + M * 2 + r vV + ^b + eA + (19) 



~ ~B M „ ■ B»» + \m 2 p h p ■ W - -F pu F^ 



g 2 p h p -W[A s g 2 s cj> 2 + A v g 2 v V^], 



(21) 



where fi^, B Mt/ and F pv are defined after eq.fl]). The 
parameters of the model are again the masses and the 
couplings, which are now constants, i.e., g s replaces T s , 
g v replaces T v and g p replaces T p . Non-linear a terms are 
also included. We have followed the prescription of 0], 
where the starting point was the NL3 parametrization 
and the g p coupling was adjusted for each value of the 
coupling Ai studied in such a way that for kp = 1.15 
fm _1 (not the saturation point) the symmetry ener gy i s 
25.68 MeV. In the present work we set A s = as in [37j |. 
Notice that other possibilities for this model with a — p 
and u — p couplings have already been discussed in the 
literature as in [4J, for instance. 

The mesonic equations of motion in the Thomas-Fermi 
approximation become 



V 2 cj ) = m 2 J-g sPs + ^cj ) 2 + ^ 3 (22) 



V 2 V Q = mlV - g vP + 2A v g 2 v V gib 2 , (23) 



Mn = ^k 2 Fn + M* 2 + T V V - ^b + £* (20) 

The values of kp p and kp n are obtained inverting these 
two last equations. 

Such density dependences in the coupling parameters 
do not affect the energy functional but of course affect its 
derivative such as the pressure density and the chemical 
otentials. As already discussed in the literature @, H, 
, [Hj], the rearrangement term is crucial in obtaining 
different behaviors in physical properties related to the 
chemical potentials or to their derivatives with respect to 
the density, such as spinodal regions, as compared with 
the more common NL3 or TM1 parametrizations. 



V 2 b Q = m 2 p b a - fp 3 + 2A v g 2 p b g 2 v V 2 , (24) 

V 2 A = -e Pp , (25) 
and the expression for the energy reads 

+ \ [(V0) 2 + m 2 s <j> 2 - (W ) 2 - m 2 v V 2 
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(Vb ) 2 -m 2 p b 2 Q ~(VA f 



m 2 v V - g v p + 2A v g 2 v V g 2 p b 2 = 0, (33) 



+g v Vop + y p 3 b + eA pp 



m 2 b - °fpz + 2A v g 2 p b g 2 v V 2 = 0, (34) 



7 +^ 



(26) 



All other expressions are very similar to the ones ob- 
tained from the TW model and can be read off from 
them bearing in mind that the density dependent cou- 
plings have to be replaced by the constant couplings. In 
particular the chemical potentials do not contain the re- 
arrangement term E^. 



IV. NUMERICAL RESULT VIA A 
NUCLEATION PROCESS 

At this point, eqs. (|12M15p for the TW model and eqs. 
(|22II25[) for the NLwp model have to be solved numerically 
in a self-consistent way and hence, initial and boundary 
conditions for each equation are necessary. One of the 
methods we use here is based on a prescription given in 
[38| , where these conditions are obtained from a situation 
of phase coexistence in a mean field approximation with 
classical meson fields and no electromagnetic interaction. 
The method is well explained in [38| and, as we are using 
different models here, just the main equations are written 
next. 

For the TW model, the equilibrium equations for ho- 
mogeneous matter for the fields are: 



m 2 (f> - 
m 2 v V a 



r s p s 
- T v p 



m pba - — P3 



0, 
0, 

0, 



(27) 
(28) 

(29) 



and for the energy and pressure density: 



p dp y p 2 + M 



*2 



+ 



+ 2 V ° + 2 °' 



(30) 



P 



—y 



7 4 p dp 



™Ld? + a^v 2 + ^b 2 

2 + 2 V ° + 2 b ° 



p dpyp 



■ M* 



+ ^ [m 2 4> 2 - m 2 v Vo - m 2 p bl] + g v V p + yp 3 &o 



6^ 24^ ~ v9v 9 p °- 



(35) 



and 



1 

3T 2 



p^dp m z s 2 mi 2 m 

--=-9 + —^V n + — On 



G 



-^ + A v glV 2 g 2 b 2 . 



(36) 



Based on the geometrical construction and Gibbs con- 
ditions for phase coexistence, i.e., the pressure and both 
chemical potentials are equal in both phases, we build 
the binodal section given in Fig. [TJ Notice that we have 
defined the proton fraction of the system as 



Vp 



Pp 
P ' 



(37) 



The binodal section yields the boundary conditions which 
we need. For the same pressure, two points, with differ- 
ent proton fractions are found. For each of these points, 
the meson fields and the densities are well defined and 
used as the initial and boundary conditions in eqs. \12\ 
[T5|) . which are then solved. Once the meson fields are 
obtained, all the quantities that depend on them, as the 
energy, pressure densities, chemical potentials, baryonic 
densities, etc are also computed. The solution is a droplet 
with a certain proton fraction surrounded by a gas of 
neutrons. If stable nuclei are calculated, the gas vanishes 
because the energy of the system lies below the neutron 
drip line and the finite nuclei properties are easily calcu- 
lated. This is the general method, but the results depend 
strongly on the model used because of the reasons dis- 
cussed in Section VI. 



+ P^o- 



(31) 

For the NLcj/? model, the equilibrium equations for ho- 
mogenous matter, energy density and pressure become: 



k A 
m 2 J - g sPs + -4> 2 + -<j> 3 = 0, 



(32) 



V. NUMERICAL RESULT WITHIN A 
HARMONIC OSCILLATOR BASIS 

Here a different prescription for solving the equations 
of motion and the thermodynamical quantities within the 
Thomas- Fermi approximation is used. According to [39l | , 
meson field equations of motion of the Klein-Gordon type 
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with sources can be carried out by an expansion in a com- 
plete set of basis states. The harmonic oscillator func- 
tions with orbital angular momentum equal to zero are 
then chosen. The oscillator length is given by 



Muj q ' 



(38) 



where M is the nucleon mass and u>o is the oscillator 
frequency. The meson fields and their corresponding in- 
homogeneous part can be expanded as 



A(r) 



TIB 



n=l 



A n R n0 {r), S A (r) = J2 S nRno{r 



where A(r) = 4>(r),Vo(r),bo(r) and 



N, 



R nl (r) = ^x l L l ^{ 2 (x 2 )exp(-x 2 /2), 
^0 



(39) 



(40) 



where x = r/bo is the radius measured in units of the 
oscillator length, 



N ni = y/2(n -!)!/(/ + n- 1/2)! 



(41) 



is the normalization constant and L™(a: 2 ) are the asso- 
ciated Laguerre polynomials. For the calculation of the 
meson fields I = in the expressions given below. Once 
the ansatz given by eqs.(|39|) are substituted into eqs.(JT2J- 
I14[) . a set of inhomogeneous equations is obtained: 



where 



H nn > = <W (bg 2 {2{n - 1) + 3/2) + m 2 A ) 



(42) 



Only the massive fields can be calculated with this 
method because the convergence of the Coulomb field, 
which has a long range, is very slow. The Green's func- 
tion method is then chosen to describe the electromag- 
netic interaction: 



A (r) = e / r^dr'p^G^r'), (44) 



with 



(45) 



, / 1/r for r > r 

GJr, r ) = < 1 c , 
v ' ' 1 1/r for r > r. 



VI. RESULTS 



Parity Violating Electron Scattering and the 
Neutron Radius 



We start this section by defining the asymmetry for 
polarized electron scattering of a hadronic target as 



A 



da+ / d£l — da-/dQ 
da + /dfl + da- 1 dfl ' 



(46) 



where da± / dD, is the differential cross section for initially 
polarized electrons with positive(+) and negative (— ) he- 
licities. As the electromagnetic interaction is not sensi- 
tive to the above difference, the asymmetry becomes de- 
pendent of the weak interaction between the electron and 
the target. Moreover, we know from the Standard Model 
that the neutral Z-boson couples more strongly to the 
neutron than to the proton. Those reasonings were then 
used in [2l| to first propose a clean way to determine 
the neutron distribution in nuclei. If we consider elastic 
scattering on an even-even target nucleus, the asymmetry 
can be written in the form: 



"i" Vnn' +\b g y/n(n + 1/2) + y/n'[n> + 1/2). 




(43) 



NL3 — - 

A y =0.01 

A V =0.025 

V 

TW 



0.1 0.2 0.3 0.4 



0.5 



A = 



Gq 1 



2naV2 



Pp(q) 



(47) 



In the above expression, G, a, a and 6& n are Standard 
Model coupling constants as defined in [2l| , q is the trans- 
ferred momentum by the electron to the nucleus and, 



n( P )(q)= d 3 r ja{qr)p n[p) {r), 



(48) 



Pn(p)( r ) being the neutron (proton) distribution in con- 
figuration space and jo the spherical Bessel function of 
order zero. It is then clear that a small q measurement 
of the asymmetry gives the neutron radius of the distri- 
bution once the proton radius is well known. The proton 
and neutron mean-square radius are defined as 



FIG. 1: Binodal section for the NL3, TW and NLwp 
parametrizations. 



2 / d 3 rr 2 Pi (r) 

tii — 



/ d 3 rpi(r) 



p,n. 



(49) 
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The neutron skin thickness is defined as 



R n — R P - 



(50) 



In the PREX experiment mentioned in the Introduc- 
tion, the asymmetry is expected to be measured at 
q w 0.4 /m _1 [22]]. Also, because the target is a heavy nu- 
cleus ( 208 Pb), the above results for the asymmetry should 
be reconsidered for a detailed comparison with the exper- 
iment, since they were obtained using a Plane Wave Born 
Approximation for the electron {43^ . For our present pur- 
poses, eq. (|4"T1) is sufficient to illustrate the sensitivity to 
the different model parametrizations and is used next in 
the presentation of our numerical results. 

The surface energy per unit area of the droplets in 
the small surface thickness approximation, excluding the 
electromagnetic field, reads [38| 



dd 
dr 



dVo 
dr 



dbo 
dr 



(51) 



However, as the electromagnetic interaction does not con- 
tribute to surface properties directly, we have kept the 
same definition for the surface energy. 

In Table II we show the neutron and proton radius, the 
neutron skin thickness, the binding energy and the sur- 
face energy obtained within the Thomas-Fermi approxi- 
mation and the two different numerical prescriptions de- 
scribed in the previous sections. All the results are sensi- 
tive to the numerical calculation although the analytical 
approximation is the same. When the nucleation method 
is performed, the neutron radius is systematically larger, 
what results in a thicker neutron skin. This is correlated 
with the fact that the surface energy is lower within the 
nucleation calculation than within the harmonic oscilla- 
tor method. Within the same numerical prescription, 
the neutron skin thickness is smaller with the TW model 
than with the NL3. As the coupling strength A„ increases 
in the NLwp model, the results move from the original 
NL3 to the TW results for all quantities, except the pro- 
ton radius, which oscillates a little. We have also in- 
cluded the results obtained with the HS parametrization 
plij because we have used this parametrization in order 
to compare the TF and the Dirac results for the cross sec- 
tions, as discussed in the following. As this parametriza- 
tion is known not to give as good results as the other 
parametrizations of the NLWM for finite nuclei, we do 
not comment on the results it provides. Notice that the 
experimental radius for the protons is obtained from the 
charge radius R c and it is given by R p = yj R% — 0.64 in 
fm [39| • Our results can be compared with experimental 
and other theoretical results found in the literature. The 
proton radius, which is known to better than 0.001 fm is 
better described within the TW model. This quantity is 
practically independent of the u> — p interaction strength 
in the NLup model as far as the HO numerical prescrip- 
tion is used. The neutron radius, on the other had, is 
strongly model dependent with drastic consequences in 



the neutron skin thickness calculation. The experimen- 
tal values for 9 are still very uncertain and all our re- 
sults fall inside the experimental confidence interval. We 
shall comment on possible restrictions to the neutron skin 
thickness in the next section. NL3 provides the best re- 
sults for the binding energy. 

In [25| . the results shown for the proton and neutron 
radius are respectively 5.52 and 5.85 fm, yielding a skin 
of 0.33 fm, larger than ours. Notice, however, that in 
[25j the Dirac equation was explicity solved. In (4j, the 
authors obtained a value of 0.21 fm for the neutron skin 
thickness and a binding energy of -7.89 MeV within a 
different parametrization of the NLujp model. Again in 
this case the Dirac equation was solved. 

In Fig. [2] we show the difference between neutron and 
proton densities at the Pb surface for the models dis- 
cussed in the present work with the Thomas-Fermi ap- 
proximation solved in a harmonic oscillator basis. While 
the curves deviate a little in between 6.0 and 8.0 fm, at 
the very surface they are similar, but a small discrepancy, 
reflecting the differences in the neutron skin can be seen. 

In Fig. [3] we display again the difference between neu- 
tron and proton densities within both numerical calcula- 
tions of the TW and NL3 models. These two Thomas- 
Fermi calculations should have given more similar results. 
However the nucleation method predicts a very small sur- 
face energy for the NL3 parametrization, and therefore, 
a large radius. This may be related to the choice of the 
boundary conditions and a deeper comparison between 
the two methods will be pursued. 

Next we present our results for the asymmetry given 
by eq. (|47[) as a function of the transfered momentum. 
We begin with Fig[4]which displays the results for the HS 
parametrization of the Walecka model. The curve labeled 
no structure means the case where Zp n (r) = Np p (r) and 
the other two curves are obtained within the TF approx- 
imation and the full solution of the Dirac equation in 



NL3 
TW 

NLoip ( A^O.01 ) 
NLtop ( A^O.025 ) 




0.000 



FIG. 2: Difference between neutron and proton densities ob- 
tained with the Thomas-Fermi approach solved in a harmonic 
oscillator basis for the models discussed in the present work. 
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the Hartree approximation. At the momentum trans- 
fer values of recent experimental interest (around ~ 0.4 
fm -1 ), the curves are almost identical. A careful analy- 
sis of the same results in a different scale shows us that 
the asymmetry changes 12 and 11 percent respectively 
within the Dirac and TF approximations in comparison 
with the no structure case. Since it is the measurement 
of the asymmetry in this low momentum transfer region 
that will provide the accurate result for the neutron skin 
thickness, we have restricted our calculations to the TF 
approximation, as stated in the Introduction. 

In Fig. [5K we show the asymmetry obtained with the 
NL3 model for both numerical calculations in the TF ap- 
proximation, i.e, nucleation and HO expansion methods. 
In this case, the agreement is very satisfactory even for 
larger g-values, although the small numerical discrepan- 
cies is reflected in a ~ 10 percent difference in the pre- 
dicted neutron skin thickness, as can be seen from Table 
II. Finally, in Fig. 03 our results for the NLup (using 
two different values for the lu — p coupling constant) and 
the TW models within the HO numerical prescription are 
shown. Again, at low momentum transfers, all curves co- 
incide. However, it should be noticed that even for two 
different model parametrizations which lead us to identi- 
cal neutron skin thicknesses, a measurement of the asym- 
metry in a higher g-region with a modest experimental 
precision, can distinguish between them. Also, we should 
expect that the asymmetry presents more structure in 
this high momentum transfer region if we solve the Dirac 
equation instead of using the TF approach, once the high 
q value region is much more sensitive to the central part 
of the neutron distribution, which is known to be flat in 
the TF approximation. These differences can be seen in 
Figg] 



VII. DIFFERENT EOS, DIFFERENT NEUTRON 
SKINS 



For the sake of completeness, at this point, we discuss 
some of the differences between the TW, the NLwp mod- 
els and the NL3 parametrization of the NLWM. 
From Fig. Q] one can see that the largest possible pres- 
sure for a phase coexistence in the TW model is much 
lower, and appears at a lower proton fraction than the 
NL3 model. This gives rise to a thinner crust within 
the TW model, which may imply that the more exotic 
pasta shapes will not form [5J. The NLup model goes 
on a different direction, i.e., the pressure becomes higher 
than the one obtained with the NL3 as the A„ coupling 
is turned on. 

Although the nuclear matter properties fitted to 
parametrize the models are quite similar (see Table , 
the way the EoS behaves when extrapolated to higher 
or lower densities can vary a lot from a density depen- 
dent hadron model to one of the parametrizations of the 
NLWM. Moreover, as seen from Table HI although the ef- 
fective mass at saturation density is lower with the TW 
than with the NL3, it can accommodate hyperons if an 
EoS for stellar matter is necessary, contrary to the usual 
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FIG. 3: Difference between neutron and proton densities ob- 
tained with the Thomas-Fermi approach solved with both nu- 
merical prescriptions for the TW model. 



FIG. 4: Parametrization HS, comparison Thomas-Fermi-HO 
versus Dirac-HO 
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NL3 parametrization [32], [33, \34 



TABLE I: Nuclear matter properties. 





NL3 








TW 




[25] 




[361 




[241 






A„ = 0.01 


A„ = 0.02 


A„ = 0.025 




B/A (MeV) 


16.3 


16.3 


16.3 


16.3 


16.3 


po (fm- 3 ) 


0.148 


0.148 


0.148 


0.148 


0.153 


K (MeV) 


271 


271 


271 


271 


240 


£sym. (MeV) 


37.4 


34.9 


33.1 


32.3 


32.0 


M*/M 


0.60 


0.60 


0.60 


0.60 


0.56 


L (MeV) 


118 


88 


68 


61 


55 


K aym (MeV) 


100 


-46 


-53 


-34 


-124 



for the TW model 0, H QJ, 13 . Again, for the sake of 
completeness we reproduce these results here because the 
neutron skin thickness and the neutron star EoS are re- 
lated by this quantity 0, [111,11, which is usually defined 

with 5 = —ps/p = 1 — 2y p . The 



„„ c _ 1 d'E/p 
ao t->sym — 2 ^2 



.5=0 



symmetry energy can be analytically rewritten as 



k 2 

F — F 
6e F 



for the TW model and as 



kp 
6e F 



r 2 

— / 

8m 2 



8m; 2 ' 



(52) 



(53) 



with the effective p-meson mass defined as [H 



Another quantity of interest in asymmetric nuclear 
matter is the nuclear bulk symmetry energy, shown in 
Table[T] for the saturation point. The differences in the 
symmetry energy at densities larger than the nuclear sat- 
uration density is still not well established, but has al- 
ready been extensively discussed in the literature even 
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for the NLwp model. In both cases 

k Fp = k F (i + 5) 1 / 3 , k Fn = Mi - <5) 1/3 , 



with k F = (1.5-7T 2 p) 1 / 3 and e F — y kp + M* 2 . In equa- 
tions (|52|) and ((53)) the second term dominates at large 
densities. It is seen that the non-linear p — uj terms 
introduce a non-linear density behavior in the symme- 
try energy of the NLWM parametrizations such as NL3 
and TM1. In TW the non- linear density behavior en- 
ters through the density dependent coupling parameters. 
These non-linear density behavior is important because 
the linear behavior of NL3 and TM1 parametrizations 
predicts too high symmetry energy at densities of impor- 
tance for neutron star matter which has direct influence 
on the proton fraction dependence with density. From 
Fig. [6l it is easily seen that the symmetry energy ob- 
tained with the TW model behaves in a very different 
way, as compared with NL3. In [4j a relation between 
the symmetry energy and the nuclear binding energy is 
discussed : the harder the EoS, the more the symmetry 
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FIG. 5: Asymmetry obtained with a) NL3 with both numer- 
ical prescriptions and b)parametrizations NLwp and TW 



FIG. 6: Symmetry energy for the NL3, TW and NLojp models. 
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energy rises with density. The density dependence dis- 
cussed in [1] is of the type introduced in [H, [3(| through 
the inclusion of a a — p and/or oj — p couplings and then, 
similar with the NLwp model discussed here. One can 
observe that as the strength of the coupling increases, 
the symmetry energy gets closer to the TW curve. In 
fact, in 8] it was shown that once this kind of coupling 
is introduced with a reasonable strength, the symmetry 
energy at low densities tends to behave as the TW model. 

The symmetry energy can be expanded around the nu- 
clear saturation density and reads 



^sym\P) — ^sym(Po) "h 



L ( P~ Po\ , K sym ( p - po 



Po 



18 



Po 



(54) 

where L and K sym are respectively the slope and the 
curvature of the nuclear symmetry energy at po and they 
are calculated from 



3po 



dp' 



„ 2 ® £sym{p) 
9P ° d*p 



(55) 

These two quantities can provide important information 
on the symmetry energy at both high and low densi- 
ties because they characterize the density dependence of 
the energy symmetry. In a recent work [49(, the authors 
found a correlation between the slope of the symmetry 
energy and the neutron skin thickness. In their work 21 
sets of the non-relativistic Skyrme potential were inves- 
tigated and only 4 of them were shown to have L values 
consistent with the values extracted from experimental 
isospin diffusion data from heavy ion collisions. In fact, 
the extracted value was L — 88 ± 25 MeV [H3|, which 
gives a very strong constraint on the density dependence 
of the nuclear symmetry energy and consequently on the 
EoS as well. A detailed analysis of Table I shows that, if 
this constraint is to be taken seriously, neither the NL3 
nor the TW model satisfy it. Nevertheless, the NLojp 
slope interpolates beautifully between the NL3 and TW 
slope values. Once again it is seen that the increase in 
A„ approximates the NL3 model values for the slope and 
energy symmetry to the TW values. Moreover, we have 
also tried to find a correlation between the 9 values shown 
in Table II and L values displayed in Table I. We found 
that, as far as some numerical imprecision are consid- 
ered, larger values of L correspond to larger values of the 
neutron skin, as seen in Fig. [7] 

Let's now go back to the problem of solving the dif- 
ferential equations within the nucleation numerical pre- 
scription. As we need boundary conditions arising from 
the liquid-gas phase coexistence in order to solve eqs. 
p2][T5)) for the TW model and eqs. (|22]|25|) for the NLwp 
model, the binodal sections are essential and the spinodal 
sections, which separate the regions of stable to unsta- 
ble matter are also of interest. If we had displayed the 
binodals in a p p versus p n plot, as it is done with the 
spinodals in Fig [51 we could see that the spinodals sur- 
faces lie inside the binodal sections and share the critical 
point corresponding to the highest pressure. 



In Fig. [8J the spinodals for the three different mod- 
els discussed in this work are shown. Once again, some 
of these results can also be found in the recent litera- 
ture @, H[ , but we include them here to make a direct 
link with the binodals. The instability of the ANM sys- 
tem is essentially determined by density fluctuations in 
the isoscalar channel. Although the spinodals are, by 
themselves, not relevant in calculations performed at the 
thermodynamical equilibrium, the isospin channel is very 
sensitive to the instabilities occurring below the nuclear 
saturation density. The spinodal is determined by the 
values of pressure, proton fraction and density for which 
the determinant of 



d 2 T 
dpidpj 



(56) 



where T is the free energy density, goes to zero. A de- 
tailed analysis of this quantity can be found in 0, . 
From Fig. [5J it is seen that the instability region in the 
p p / p n plane, defined by the inner section of the spinodal 
curve is larger for the TW than for the NL3 model. The 
size of the instability region depends on the derivative 
of the chemical potentials with respect to the neutron 
and proton densities. At low densities different models 
exhibit different behaviors. 

The presence of the rearrangement term in the TW 
model also plays a decisive role. Even though a rela- 
tively large compensation exists between scalar and vec- 
tor mesons in the isoscalar channels within the rearrange- 
ment term at low densities, the spinodal region is defined 
by the derivative of the chemical potential and therefore 
of the rearrangement term. 

Next we examine the spinodals obtained with differ- 
ent coupling strengths for the NLujp model. As seen in 
Fig. [5J there is almost no difference between the different 
curves. They all fall around the original NL3 curve but 
once again, they tend to the TW curve as the coupling 
strength increases. However, contrary to the TW model, 
it was shown in 0] that the direction of the instability in 
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FIG. 7: Correlation between the neutron skin 9 and the slope 
of the symmetry energy L. 



11 



NL W(3 increases distillation as the density increases, and 
the larger the coupling A v the larger the effect. 

Finally, to end this section, let's make our points clear: 
we have used a simple mean field theory approach to ob- 
tain the boundary conditions for the equations of motion 
of the meson fields in the nucleation prescription. These 
boundary conditions depend on the model used and are 
intrinsically related with the liquid-gas phase transition 
which, in turn, can be well understood by studying the 
coexistence surfaces of the corresponding models. On 
the other hand, the neutron skin thickness shows a lin- 
ear correlation with the slope of the symmetry energy, 
as already pointed out in [491 ] for non-relativistic mod- 
els. Based on the different behaviors found with density 
dependent hadronic models and the NLWM, an obvious 
consequence is the fact that the neutron skin thickness 
depends on the choice of the model. 

VIII. CONCLUSIONS 

We have calculated the 208 Pb neutron skin thickness 
with two different density dependent hadronic models, 
the TW and the NLwp model, and one of the most used 
parametrizations of the NLWM, the NL3. The calcu- 
lations were done within the Thomas-Fermi approxima- 
tion, which gives quite accurate results for the asymme- 
try in the momentum transfer range of interest for the 
calculation of neutron skins. In implementing the nu- 
merical results two different prescriptions were used: the 
first one based on the nucleation process and the second 
one based on the harmonic oscillator basis method. We 
have seen that when the nucleation method is performed, 
the neutron radius is systematically larger, what results 



in a thicker neutron skin. This is a consequence of the 
fact that the surface energy is lower within the nucleation 
calculation than within the harmonic oscillator method. 
Within the same numerical prescription, the neutron skin 
thickness is smaller with the TW model than with the 
NL3. As the coupling strength A„ increases in the NLwp 
model, the neutron skin thickness moves from the orig- 
inal NL3 towards the TW results. We have also found 
that although the neutron skin thickness is model depen- 
dent, the asymmetry at low momentum transfers (below 
0.5 fm -1 ) is very similar for all models and all numerical 
prescriptions. As q increases, the asymmetry also be- 
comes model dependent. The density profiles obtained 
from the solution of the Dirac equation exhibits oscil- 
lations near the center of the nucleus, behavior which is 
not reproduced within the Thomas- Fermi approximation. 
This fact shows up in the asymmetry at large momentum 
transfers and therefore all the calculations should be re- 
produced by solving the Dirac equation. This calculation 
is already under investigation. 



It is worth mentioning that the neutron skin thickness 
has shown to give hints on the equations of state that 
are suitable to describe neutron stars. Moreover, in (4|| 
a correlation between the slope of the symmetry energy 
and the neutron skin thickness was found for Skyrme- 
type models. We have observed that this correlation was 
also present in the density dependent models we have 
studied in the present work. 
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FIG. 8: Spinodal section in terms of p p versus p n for the NL3, 
TW and NLujp models. 
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TABLE II: 208 Pb properties 
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